#include <arm_neon.h>
#include <pthread.h>
#include <stdlib.h>
#include <iostream>
#include <sys/time.h>
#include <semaphore.h>
#include <unistd.h>
#include <omp.h>
using namespace std;
#define NUM_THREADS 7
//线程数据结构定义
typedef struct {

    int k; //消去的轮次
    int t_id; // 线程 id
    float** a;
    int n;
    int all;
}threadParam_t;
//信号量定义
sem_t sem_main;
sem_t sem_workerstart[NUM_THREADS]; // 每个线程有自己专属的信号量
sem_t sem_workerend[NUM_THREADS];
void Serial(float** a, int n)
{
    int i, j, k;
    for (k = 0; k < n; k++)
    {
        for (j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        for (i = k + 1; i < n; i++)
        {
            for (j = k + 1; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0;
        }
    }
}
void* State_threadFunc_Serial(void* param)
{
    threadParam_t* p = (threadParam_t*)param;
    int t_id = p->t_id;
    float** a = p->a;
    int n = p->n;
    for (int k = 0; k < n; ++k)
    {
        sem_wait(&sem_workerstart[t_id]); // 阻塞，等待主线完成除法操作
        //循环划分任务
        for (int i = k + 1 + t_id; i < n; i += NUM_THREADS)
        {
            //消去
            for (int j = k + 1; j < n; ++j)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0.0;
        }
        sem_post(&sem_main); // 唤醒主线程
        sem_wait(&sem_workerend[t_id]); //阻塞，等待主线程唤醒进入下一轮
    }
    pthread_exit(NULL);
}
void State_Pthread_Serial(float** a, int n)
{
    //初始化信号量
    sem_init(&sem_main, 0, 0);
    for (int i = 0; i < NUM_THREADS; ++i)
    {
        sem_init(&sem_workerstart[i], 0, 0);
        sem_init(&sem_workerend[i], 0, 0);
    }
    //创建线程
    pthread_t* handles = (pthread_t*)malloc(NUM_THREADS * sizeof(pthread_t));// 创建对应的 Handle
    threadParam_t* param = (threadParam_t*)malloc(NUM_THREADS * sizeof(threadParam_t));// 创建对应的线程数据结构
    for (int t_id = 0; t_id < NUM_THREADS; t_id++)
    {
        param[t_id].t_id = t_id;
        param[t_id].k = 0;
        param[t_id].a = a;
        param[t_id].n = n;
        param[t_id].all = 0;
        pthread_create(handles + t_id, 0, State_threadFunc_Serial, &param[t_id]);
    }
    for (int k = 0; k < n; ++k)
    {
        //主线程做除法操作
        for (int j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        //开始唤醒工作线程
        for (int t_id = 0; t_id < NUM_THREADS; ++t_id)
            sem_post(&sem_workerstart[t_id]);
        //主线程睡眠（等待所有的工作线程完成此轮消去任务）
        for (int t_id = 0; t_id < NUM_THREADS; ++t_id)
            sem_wait(&sem_main);
        // 主线程再次唤醒工作线程进入下一轮次的消去任务
        for (int t_id = 0; t_id < NUM_THREADS; ++t_id)
            sem_post(&sem_workerend[t_id]);
    }
    for (int t_id = 0; t_id < NUM_THREADS; t_id++)
        pthread_join(handles[t_id], 0);
    sem_destroy(&sem_main);
    for (int i = 0;i < NUM_THREADS;i++)
    {
        sem_destroy(&sem_workerend[i]);
        sem_destroy(&sem_workerstart[i]);
    }
}
void* Dynamic_threadFunc_Serial(void* param)
{
    threadParam_t* p = (threadParam_t*)param;
    int k = p->k; //消去的轮次
    int t_id = p->t_id; //线程编号
    float** a = p->a;
    int n = p->n;
    int all = p->all;
    int i = k + t_id * all + 1; //获取自己的计算任务
    for (; i < k + t_id * all + 1 + all; i++)
    {
        for (int j = k + 1; j < n; ++j)
            a[i][j] = a[i][j] - a[i][k] * a[k][j];
        a[i][k] = 0;
    }
    pthread_exit(NULL);
}
void Dynamic_Pthread_Serial(float** a, int n)
{
    for (int k = 0; k < n; ++k)
    {
        //主线程做除法操作
        for (int j = k + 1; j < n; j++)
            a[k][j] = a[k][j] / a[k][k];
        a[k][k] = 1.0;
        //创建工作线程，进行消去操作
        int worker_count = 7; //工作线程数量
        pthread_t* handles = (pthread_t*)malloc(worker_count * sizeof(pthread_t));// 创建对应的 Handle
        threadParam_t* param = (threadParam_t*)malloc(worker_count * sizeof(threadParam_t));// 创建对应的线程数据结构
        //分配任务
        int t = (n - k - 1) / worker_count;
        for (int t_id = 0; t_id < worker_count; t_id++)
        {
            param[t_id].k = k;
            param[t_id].t_id = t_id;
            param[t_id].a = a;
            param[t_id].n = n;
            param[t_id].all = t;
        }
        //创建线程
        for (int t_id = 0; t_id < worker_count; t_id++)
            pthread_create(&handles[t_id], 0, Dynamic_threadFunc_Serial, &param[t_id]);
        //主线程挂起等待所有的工作线程完成此轮消去工作
        for (int t_id = 0; t_id < worker_count; t_id++)
            pthread_join(handles[t_id], 0);
        for (int i = k + t * worker_count + 1; i < n; i++)
        {
            for (int j = k + 1; j < n; j++)
                a[i][j] = a[i][j] - a[i][k] * a[k][j];
            a[i][k] = 0;
        }
    }
}
void OpenMP_Static(float** a, int n)
{
#pragma omp parallel num_threads(8) shared(a,n)
    {
        for (int k = 0; k < n; k++)
        {
#pragma omp for schedule(static)
            for (int j = k + 1; j < n; j++)
                a[k][j] = a[k][j] / a[k][k];
#pragma omp barrier
            ;
#pragma omp single
            a[k][k] = 1.0;
#pragma omp barrier
            ;
#pragma omp for collapse(2) schedule(static)
            for (int i = k + 1; i < n; i++)
                for (int j = k + 1; j < n; j++)
                    a[i][j] = a[i][j] - a[i][k] * a[k][j];
#pragma omp barrier
            ;
#pragma omp for schedule(static)
            for (int i = k + 1; i < n; i++)
                a[i][k] = 0;
        }
    }
}
void OpenMP_Dynamic(float** a, int n)
{
#pragma omp parallel num_threads(8) shared(a,n)
    {
        for (int k = 0; k < n; k++)
        {
#pragma omp for schedule(guided)
            for (int j = k + 1; j < n; j++)
                a[k][j] = a[k][j] / a[k][k];
#pragma omp barrier
            ;
#pragma omp critical
            a[k][k] = 1.0;
#pragma omp barrier
            ;
#pragma omp for collapse(2) schedule(guided)
            for (int i = k + 1; i < n; i++)
                for (int j = k + 1; j < n; j++)
                    a[i][j] = a[i][j] - a[i][k] * a[k][j];
#pragma omp barrier
            ;
#pragma omp for schedule(dynamic)
            for (int i = k + 1; i < n; i++)
                a[i][k] = 0;
        }
    }
}
void OpenMP_Guided(float** a, int n)
{
#pragma omp parallel num_threads(8) shared(a,n)
    {
        for (int k = 0; k < n; k++)
        {
#pragma omp for schedule(guided)
            for (int j = k + 1; j < n; j++)
                a[k][j] = a[k][j] / a[k][k];
#pragma omp barrier
            ;
#pragma omp critical
            a[k][k] = 1.0;
#pragma omp barrier
            ;
#pragma omp for collapse(2) schedule(guided)
            for (int i = k + 1; i < n; i++)
                for (int j = k + 1; j < n; j++)
                    a[i][j] = a[i][j] - a[i][k] * a[k][j];
#pragma omp barrier
            ;
#pragma omp for schedule(guided)
            for (int i = k + 1; i < n; i++)
                a[i][k] = 0;
        }
    }
}
void m_rest(float** A, int n)
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < i; j++)
            A[i][j] = 0;
        A[i][i] = 1.0;
        for (int j = i + 1; j < n; j++)
        {
            float x = rand();
            while (x == 0) x = rand();
            A[i][j] = x;
        }
    }
    for (int k = 0; k < n; k++)
        for (int i = k + 1; i < n; i++)
            for (int j = 0; j < i; j++)
                A[i][j] += A[k][j];
}


int main()
{
    srand((unsigned)time(NULL));
    int n = 10;
    while (n < 1600)
    {
        float** a = (float**)aligned_alloc(32, n * sizeof(float*));
        for (int i = 0; i < n; i++)
            a[i] = (float*)aligned_alloc(32, n * sizeof(float));
        m_rest(a, n);

        cout << n << " ";
        struct timeval tstart, tend;
        double timeUsed = 0.0;
        float** t = (float**)aligned_alloc(32, n * sizeof(float*));
        for (int i = 0; i < n; i++)
            t[i] = (float*)aligned_alloc(32, n * sizeof(float));

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        Serial(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        double time1 = timeUsed;
        cout << timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        State_Pthread_Serial(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << time1/timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        Dynamic_Pthread_Serial(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << time1 / timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        OpenMP_Static(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << time1 / timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        OpenMP_Dynamic(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << time1 / timeUsed << " ";

        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                t[i][j] = a[i][j];
        gettimeofday(&tstart, NULL);
        OpenMP_Guided(t, n);
        gettimeofday(&tend, NULL);
        timeUsed = 1000000 * (tend.tv_sec - tstart.tv_sec) + tend.tv_usec - tstart.tv_usec;
        timeUsed /= 1000;
        cout << time1 / timeUsed << " ";



        cout << endl;
        for (int i = 0; i < n; i++)
        {
            delete[] t[i];
            delete[] a[i];
        }
        delete[] t;
        delete[] a;
        if (n >= 800) n += 200;
        if (n >= 400 && n < 800) n += 100;
        if (n >= 200 && n < 400) n += 50;
        if (n < 200) n += 10;
    }
    return 0;

}
